projekt2

library(dplyr)

## 
## Dołączanie pakietu: 'dplyr'

## Następujące obiekty zostały zakryte z 'package:stats':
## 
##     filter, lag

## Następujące obiekty zostały zakryte z 'package:base':
## 
##     intersect, setdiff, setequal, union

library(Hmisc)

## Ładowanie wymaganego pakietu: lattice

## Ładowanie wymaganego pakietu: survival

## Ładowanie wymaganego pakietu: Formula

## Ładowanie wymaganego pakietu: ggplot2

## 
## Dołączanie pakietu: 'Hmisc'

## Następujące obiekty zostały zakryte z 'package:dplyr':
## 
##     src, summarize

## Następujące obiekty zostały zakryte z 'package:base':
## 
##     format.pval, units

library(fitdistrplus)

## Ładowanie wymaganego pakietu: MASS

## 
## Dołączanie pakietu: 'MASS'

## Następujący obiekt został zakryty z 'package:dplyr':
## 
##     select

library(logspline)
library(Sim.DiffProc)

## Package 'Sim.DiffProc', version 4.8
## browseVignettes('Sim.DiffProc') for more informations.

library(gap)

## Ładowanie wymaganego pakietu: gap.datasets

## gap version 1.2.3-6

library(DescTools)

## 
## Dołączanie pakietu: 'DescTools'

## Następujące obiekty zostały zakryte z 'package:Sim.DiffProc':
## 
##     Median, Mode

## Następujące obiekty zostały zakryte z 'package:Hmisc':
## 
##     %nin%, Label, Mean, Quantile

library(gamlss)

## Ładowanie wymaganego pakietu: splines

## Ładowanie wymaganego pakietu: gamlss.data

## 
## Dołączanie pakietu: 'gamlss.data'

## Następujący obiekt został zakryty z 'package:datasets':
## 
##     sleep

## Ładowanie wymaganego pakietu: gamlss.dist

## 
## Dołączanie pakietu: 'gamlss.dist'

## Następujący obiekt został zakryty z 'package:Sim.DiffProc':
## 
##     BB

## Ładowanie wymaganego pakietu: nlme

## 
## Dołączanie pakietu: 'nlme'

## Następujący obiekt został zakryty z 'package:dplyr':
## 
##     collapse

## Ładowanie wymaganego pakietu: parallel

##  **********   GAMLSS Version 5.4-3  **********

## For more on GAMLSS look at https://www.gamlss.com/

## Type gamlssNews() to see new features/changes/bug fixes.

## 
## Dołączanie pakietu: 'gamlss'

## Następujący obiekt został zakryty z 'package:gap':
## 
##     cs

library(gamlss.dist)
library(gamlss.add)

## Ładowanie wymaganego pakietu: mgcv

## This is mgcv 1.8-40. For overview type 'help("mgcv-package")'.

## Ładowanie wymaganego pakietu: nnet

## 
## Dołączanie pakietu: 'nnet'

## Następujący obiekt został zakryty z 'package:mgcv':
## 
##     multinom

## Ładowanie wymaganego pakietu: rpart

library(hash)

## hash-2.2.6.2 provided by Decision Patterns

library(reshape)

## 
## Dołączanie pakietu: 'reshape'

## Następujący obiekt został zakryty z 'package:dplyr':
## 
##     rename

library(plotly)

## 
## Dołączanie pakietu: 'plotly'

## Następujący obiekt został zakryty z 'package:reshape':
## 
##     rename

## Następujący obiekt został zakryty z 'package:MASS':
## 
##     select

## Następujący obiekt został zakryty z 'package:Hmisc':
## 
##     subplot

## Następujący obiekt został zakryty z 'package:ggplot2':
## 
##     last_plot

## Następujący obiekt został zakryty z 'package:stats':
## 
##     filter

## Następujący obiekt został zakryty z 'package:graphics':
## 
##     layout

library(glmnet)

## Ładowanie wymaganego pakietu: Matrix

## 
## Dołączanie pakietu: 'Matrix'

## Następujący obiekt został zakryty z 'package:reshape':
## 
##     expand

## Loaded glmnet 4.1-4

library(caret)

## 
## Dołączanie pakietu: 'caret'

## Następujący obiekt został zakryty z 'package:gamlss':
## 
##     calibration

## Następujące obiekty zostały zakryte z 'package:DescTools':
## 
##     MAE, RMSE

## Następujący obiekt został zakryty z 'package:survival':
## 
##     cluster

library(randomForest)

## randomForest 4.7-1.1

## Type rfNews() to see new features/changes/bug fixes.

## 
## Dołączanie pakietu: 'randomForest'

## Następujący obiekt został zakryty z 'package:ggplot2':
## 
##     margin

## Następujący obiekt został zakryty z 'package:dplyr':
## 
##     combine

library(sail)
library(splitTools)
library(groupdata2)

## 
## Dołączanie pakietu: 'groupdata2'

## Następujący obiekt został zakryty z 'package:splitTools':
## 
##     partition

Zadanie 1

Wczytujemy dane

X_train <- read.csv(file = './data/X_train.csv')
y_train <- read.csv(file = './data/y_train.csv')
X_test <- read.csv(file = './data/X_test.csv')

dim(X_train)

## [1] 3794 9000

dim(X_test)

## [1]  670 9000

dim(y_train)

## [1] 3794    1

Dane treningowe zawierają 3794 obserwacji, z kolei dane testowe 670. Zmiennych jest 9000.

Przyjrzyjmy się typom danych.

length(select_if(X_train,is.numeric))

## [1] 9000

length(select_if(X_test,is.numeric))

## [1] 9000

str(y_train)

## 'data.frame':    3794 obs. of  1 variable:
##  $ CD36: num  1.29 2.26 1.78 2 0 ...

Z opisu danych jasno wynika, że powinniśmy się spodziewać danych w pełni numerycznych. Faktycznie tak jest, każda z kolumn jest typu numerycznego.

Sprawdźmy czy dane są kompletne.

sum(is.na(X_train))

## [1] 0

sum(is.na(X_test))

## [1] 0

sum(is.na(y_train))

## [1] 0

Dane są kompletne.

Policzmy najpierw kilka podstawowych statystyk i pokażmy histogram

summary(y_train)

##       CD36       
##  Min.   :0.0000  
##  1st Qu.:0.2874  
##  Median :1.3865  
##  Mean   :1.1229  
##  3rd Qu.:1.7828  
##  Max.   :2.8647

y_train.vector <- unlist(as.vector(y_train))
hist(y_train.vector,breaks = 45)

Na ten moment ciężko stwierdzić, jaki to rozkład.

Użyjemy teraz funkcji descdist, żeby znaleźć kandydatów na poszukiwany rozkład.

descdist(y_train.vector, discrete = FALSE, boot=500)

## summary statistics
## ------
## min:  0   max:  2.864712 
## median:  1.386459 
## mean:  1.122854 
## estimated sd:  0.7749966 
## estimated skewness:  -0.3154322 
## estimated kurtosis:  1.591972

Z wykresu wynika, że kandydatami na rozkład jest rozkład jestostajny lub rozkład beta.

Sprawdźmy jak wygląda rozkład epmiryczny w porównaniu do rozkładu jednostajnego.

fit.unif <- fitdist(y_train.vector, "unif")
fit.unif

## Fitting of the distribution ' unif ' by maximum likelihood 
## Parameters:
##     estimate Std. Error
## min 0.000000         NA
## max 2.864712         NA

plot(fit.unif)

Naszą uwagę w szczególności powinien przykuć QQ-plot. Widzimy, że moglibyśmy z dużym przymróżeniem oka przyjąć, że wykres nawija się na linię prostą, jednak początek wykresu może nas w szczególności martwić. Również histogram nie fituje zbyt dobrze.

Sprawdźmy zatem rozkład beta.

Niestety nasz rozkład nie jest w [0,1], w przeciwieństwie do rozkładu beta. Zatem wstrzymujemy nasze rozważania nad tym kandydatem.

Weźmy jeszcze dla porównania rozkład normalny.

fit.norm <- fitdist(y_train.vector, "norm")
plot(fit.norm)

Tutaj QQ-plot wygląda koszmarnie, sytuacja na początku wykresu jest jeszcze gorsza niż w przypadku uniform. Histogram też nie fituje zbyt dobrze.

Zatem na podstawie tych wykresów póki co optowalibyśmy za rozkładem jednostajnym

Sprawdźmy teraz Kryterium informacyjne Akaikego (AIC)

fit.norm$aic

## [1] 8835.75

fit.unif$aic

## [1] 7990.126

Dla potwierdzenia, że rozkład jednostajny pasuje tutaj lepiej, widzimy, że aic jest dla niego mniejsze, niż dla rozkładu normalnego.

Przeprowadźmy teraz test Kolmogorova-Smirnova

set.seed(5)

sample.unif.distribution <- runif(3794, min=0, max=2.864712 ) # tworzymy wektor z rozkładu jednostajnego o wyestymowanych wcześniej min i max

ks.test(y_train.vector, sample.unif.distribution)

## Warning in ks.test.default(y_train.vector, sample.unif.distribution): wartość
## prawdopodobieństwa w obecności powtórzonych wartości będzie przybliżona

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  y_train.vector and sample.unif.distribution
## D = 0.22114, p-value < 2.2e-16
## alternative hypothesis: two-sided

p-value jest bardzo małe, więc odrzucamy hipotezę zerową, że rozkłady są równe.

Sprawdźmy jeszcze jak wygląda wykres ECDF vs CDFs np. dla rozkładu normalnego, jednostajnego, czy dla rozkładu exp.

fit.norm <- fitdist(y_train.vector, "norm")
fit.exp <- fitdist(y_train.vector, "exp")
fit.unif <- fitdist(y_train.vector, "unif")


cdfcomp(list(fit.norm, fit.exp, fit.unif),
 legendtext = c("norm", "exp", "unif"))

Widzimy, że żaden z rozkładów nie fituje jakoś specjalnie dobrze.

Spróbujmy innej metody. fitDist z biblioteki fitdistrplus stara się dopasować nasz rozkład do wielu znanych rozkładów, wybierając najlepszy według mle.

fit <- fitDist(y_train.vector, k = 2, type = "realplus", trace = FALSE, try.gamlss = TRUE)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |======                                                                |   9%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |=========                                                             |  13%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |===========================                                           |  39%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |==============================                                        |  43%

## Warning in MLE(ll2, start = list(eta.mu = eta.mu, eta.sigma = eta.sigma), :
## possible convergence problem: optim gave code=1 false convergence (8)

## 
  |                                                                            
  |=================================                                     |  48%

## Warning in MLE(ll2, start = list(eta.mu = eta.mu, eta.sigma = eta.sigma), :
## possible convergence problem: optim gave code=1 false convergence (8)

## 
  |                                                                            
  |=====================================                                 |  52%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |========================================                              |  57%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |===========================================                           |  61%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |==============================================                        |  65%

## Warning in MLE(ll3, start = list(eta.mu = eta.mu, eta.sigma = eta.sigma, :
## possible convergence problem: optim gave code=1 false convergence (8)

## 
  |                                                                            
  |=================================================                     |  70%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |====================================================                  |  74%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |================================================================      |  91%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |===================================================================   |  96%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## 
  |                                                                            
  |======================================================================| 100%

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

## Warning in nlminb(start = start, objective = f, control = optim.control): NA/NaN
## function evaluation

summary(fit)

## *******************************************************************
## Family:  c("EXP", "Exponential") 
## 
## Call:  gamlssML(formula = y, family = EXP) 
## 
## Fitting method: "nlminb" 
## 
## 
## Coefficient(s):
##         Estimate  Std. Error  t value   Pr(>|t|)    
## eta.mu  0.115873    0.016235  7.13727 9.5191e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Degrees of Freedom for the fit: 1 Residual Deg. of Freedom   3793 
## Global Deviance:     8467.25 
##             AIC:     8469.25 
##             SBC:     8475.49

Najlepszym otrzymanym rozkładem okazał się exp(0.115873), ponownie spróbujmy testu Kolmogorova-Smirnova

sample.exp.distribution <- rexp(3794, rate = 0.115873) # tworzymy wektor z rozkładu wykładniczego z wyestymwoanym parametrem

ks.test(y_train.vector, sample.exp.distribution)

## Warning in ks.test.default(y_train.vector, sample.exp.distribution): wartość
## prawdopodobieństwa w obecności powtórzonych wartości będzie przybliżona

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  y_train.vector and sample.exp.distribution
## D = 0.76173, p-value < 2.2e-16
## alternative hypothesis: two-sided

Tutaj również p-value bardzo małe, więc odrzucamy hipotezę zerową, że rozkłady są równe.

Na podstawie powyższych rozważań nie udało nam się dopasować rozkładu empirycznego do żadnego znanego nam rozkładu.

Wybierzemy najpierw 250 najbardziej skorelowanych zmiennych ze zmienną objaśnianą. Napiszmy funkcję…

get_most_correlated <- function(x,y,n) {
  # funkcja zwracająca nazwy n najbardziej skorelowanych zmiennych ze zmienną objaśnianą
  # input: x - DataFrame odpowiadający zmiennym objaśniającym, y - dataframe odpowiadający zmiennej objaśnianej,
  #        n - zmienna typu integer ile skorelowanych zmiennych chcemy zwrócić
  # output: lista nazw zmiennych
  
  h <- list() # w tej liście przechowamy informację o korelacjach
  
  for (i in 1:ncol(x))
  {
      h[[colnames(x[i])]] <- abs(cor.test(unlist(as.vector(x[i])), unlist(as.vector(y)), method="pearson")$estimate)
      # dla każdej zmiennej wprowadzamy wartość jej korelacji ze zmienną objaśnianą do listy typu key-value
  }
    
  h <- h[order(unlist(h), decreasing=TRUE)] # sortujemy listę malejąco według wartości
  
  return(names(h[1:n])) # zwracamy n pierwszych wartości

}

most_correlated <- get_most_correlated(X_train,y_train,250) # zgarniamy 250 najardziej skorelowanych wartości

Policzmy korelacje między parami oraz wyświetlmy heatmapę. Myślę, że najlepszą heatmapę w naszym przypadku uzyskamy korzystając z biblioteki plotly.

corr <- cor(X_train[,c(most_correlated)]) # zgarniamy korelacje między sobą pomiędzy uzyskanymi zmiennymi
corr <- round(corr, 2) # wyniki zaokrąglamy

plot_ly(z = corr, x = most_correlated, y = most_correlated, type = "heatmap")   

Zadanie 2

1. Elastic Net pojawił się po raz pierwszy w wyniku krytyki lassa, którego dobór zmiennych może być zbyt zależny od danych, a przez to niestabilny. Regresja Elastic Net łączy w sobie cechy Ridge oraz Lasso, poprzez zastosowanie dwóch typów norm, których istotność względem siebie kontroluje hiperparametr α. Elastic Net najlepiej nadaje się do modelowania danych z dużą liczbą wysoce skorelowanych predyktorów. Warto zaznaczyć, że pomimo ogromnej ilości zmiennych nie zajmiemy się future selection, bo to LASSO, czy elastic net dokonują feature selection, nakładając na odpowiednie cechy współczynnik = 0.

elastic net 1

Zadaniem elastic net jest zminimalizowanie następującej funkcji:

Uzyskując przy tym estymator wektora parametrów β.

Hiperparametry: - λ ( parametr regularyzacji ), to stała mnożąca warunki kary. λ = 0 jest odpowiednikiem zwykłej metody najmniejszych kwadratów, rozwiązywanej przez regresję liniową. Z kolei λ = 1 jest używana do otrzymania pełnej kary. Gdy λ zwiększa się w nieskończoność, to efekt regularyzacji jest wzmocniony i jedynym celem staje się utrzymanie małych współczynników β, zamiast minimalizowania funkcji straty. Powszechne są bardzo małe wartości lambady, takie jak 1e-3 lub mniejsze. - α ( parametr mieszania ), należy do przedziału [0,1]. Dla α = 0 otrzymujemy regresję grzbietową, a dla α = 1 otrzymujemy LASSO.

Tworzymy grid używając funkcji z biblioteki caret.

Uwaga: naszym kryterium błędu będzie pierwiastek ze sredniego błędu kwadratowego.

grid <- expand.grid(alpha = c(0, 0.2, 0.4, 0.6, 0.8, 1),
 lambda = c(0.001, 0.01, 0.1, 1, 10))

grid

##    alpha lambda
## 1    0.0  1e-03
## 2    0.2  1e-03
## 3    0.4  1e-03
## 4    0.6  1e-03
## 5    0.8  1e-03
## 6    1.0  1e-03
## 7    0.0  1e-02
## 8    0.2  1e-02
## 9    0.4  1e-02
## 10   0.6  1e-02
## 11   0.8  1e-02
## 12   1.0  1e-02
## 13   0.0  1e-01
## 14   0.2  1e-01
## 15   0.4  1e-01
## 16   0.6  1e-01
## 17   0.8  1e-01
## 18   1.0  1e-01
## 19   0.0  1e+00
## 20   0.2  1e+00
## 21   0.4  1e+00
## 22   0.6  1e+00
## 23   0.8  1e+00
## 24   1.0  1e+00
## 25   0.0  1e+01
## 26   0.2  1e+01
## 27   0.4  1e+01
## 28   0.6  1e+01
## 29   0.8  1e+01
## 30   1.0  1e+01

Teraz zajmiemy się wyborem liczby podzbiorów walidacji krzyżowej. Standardowym wyborem ilości podzbiorów jest k = 10 lub k = 5. Jednak k=10 wybiera się z reguły wtedy gdy ilość danych jest duża, a nie jest to nasz przypadek. Zatem raczej optowalibyśmy za k = 5, czy też k = 4. W ogólności wybór nie jest ściśle określony, bo trudno oszacować jak dobrze nasz podzbiór tzw. fold, reprezentuje cały zbiór danych. k = 5 oznaczałoby, że każdy podzbiór składałby się z około 3790 * 20% = 758 obserwacji. Taka ilość danych wydaje się sensowna dlatego zdecydujemy się na k = 5.

Gdyby nasz dataset treningowy obejmowal np. 100000 obserwacji, wtedy dla k=10 pojedynczy fold odnosiłby się do 10000 obserwacji, co w ogólności wystarcza na wykonanie rzetelnego testu. Dodatkowo nie skupiamy się tu na aspekcie obliczeń numerycznych, bo mamy dostatecznie dobre komputery żeby różnica dla naszych danych pomiędzy k = 5 i k = 10 nie miała większego znaczenia jeśli chodzi o czas obliczeń.

Najpierw przeskalujmy dane.

set.seed(1)

elnet.models <- list() # tu dodamy kolejne modele elastic net

X_train_scaled = scale(X_train)
X_test_scaled = scale(X_test, center=attr(X_train_scaled, "scaled:center"), 
                              scale=attr(X_train_scaled, "scaled:scale"))

X_train_scaled <- data.frame(X_train_scaled)
X_test_scaled <- data.frame(X_test_scaled)


folds <- split(X_train_scaled, sample(1:5, nrow(X_train_scaled), replace=T)) # dzielimy dataset na 5 podzbiorów

ind <- list() # indeksy kolejnych podzbiorów

for (j in 1:5){
  
  ind <- append(ind,list(rownames(data.frame(folds[j])))) # zgarniamy indeksy kolejnych podzbiorów
  
}

for (j in 1:5){
  print('.')
  elnet <- train(x=X_train_scaled[unlist(ind[j]),], y= y_train[as.numeric(unlist(ind[j])),], method = "glmnet",  tuneGrid=as.data.frame(grid))
  print('..')
  elnet.models <- list(elnet.models, elnet)
}

## [1] "."

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.

## [1] ".."
## [1] "."

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.

## [1] ".."
## [1] "."

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.

## [1] ".."
## [1] "."

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.

## [1] ".."
## [1] "."

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.

## [1] ".."

#elnet.models

Liczymy foldsy do walidacji krzyżowej.

flds <- createFolds(y=unlist(y_train), k = 5, list = TRUE, returnTrain = FALSE)

train.control <- trainControl( 
    index=flds
    , verboseIter = T 
    , returnData = T 
    , savePredictions = T 
    ) 

set.seed(1)

X_train_scaled = scale(X_train)
X_test_scaled = scale(X_test, center=attr(X_train_scaled, "scaled:center"), 
                              scale=attr(X_train_scaled, "scaled:scale"))

X_train_scaled <- data.frame(X_train_scaled)
X_test_scaled <- data.frame(X_test_scaled)

#cv_5 <- trainControl(method = "cv", number = 5)

y_train <- as.numeric(unlist(y_train))
elnet <- train(x=X_train_scaled, y=y_train, method = "glmnet",trControl = train.control, tuneGrid=as.data.frame(grid))

## + Fold1: alpha=0.0, lambda=10 
## - Fold1: alpha=0.0, lambda=10 
## + Fold1: alpha=0.2, lambda=10 
## - Fold1: alpha=0.2, lambda=10 
## + Fold1: alpha=0.4, lambda=10 
## - Fold1: alpha=0.4, lambda=10 
## + Fold1: alpha=0.6, lambda=10 
## - Fold1: alpha=0.6, lambda=10 
## + Fold1: alpha=0.8, lambda=10 
## - Fold1: alpha=0.8, lambda=10 
## + Fold1: alpha=1.0, lambda=10 
## - Fold1: alpha=1.0, lambda=10 
## + Fold2: alpha=0.0, lambda=10 
## - Fold2: alpha=0.0, lambda=10 
## + Fold2: alpha=0.2, lambda=10 
## - Fold2: alpha=0.2, lambda=10 
## + Fold2: alpha=0.4, lambda=10 
## - Fold2: alpha=0.4, lambda=10 
## + Fold2: alpha=0.6, lambda=10 
## - Fold2: alpha=0.6, lambda=10 
## + Fold2: alpha=0.8, lambda=10 
## - Fold2: alpha=0.8, lambda=10 
## + Fold2: alpha=1.0, lambda=10 
## - Fold2: alpha=1.0, lambda=10 
## + Fold3: alpha=0.0, lambda=10 
## - Fold3: alpha=0.0, lambda=10 
## + Fold3: alpha=0.2, lambda=10 
## - Fold3: alpha=0.2, lambda=10 
## + Fold3: alpha=0.4, lambda=10 
## - Fold3: alpha=0.4, lambda=10 
## + Fold3: alpha=0.6, lambda=10 
## - Fold3: alpha=0.6, lambda=10 
## + Fold3: alpha=0.8, lambda=10 
## - Fold3: alpha=0.8, lambda=10 
## + Fold3: alpha=1.0, lambda=10 
## - Fold3: alpha=1.0, lambda=10 
## + Fold4: alpha=0.0, lambda=10 
## - Fold4: alpha=0.0, lambda=10 
## + Fold4: alpha=0.2, lambda=10 
## - Fold4: alpha=0.2, lambda=10 
## + Fold4: alpha=0.4, lambda=10 
## - Fold4: alpha=0.4, lambda=10 
## + Fold4: alpha=0.6, lambda=10 
## - Fold4: alpha=0.6, lambda=10 
## + Fold4: alpha=0.8, lambda=10 
## - Fold4: alpha=0.8, lambda=10 
## + Fold4: alpha=1.0, lambda=10 
## - Fold4: alpha=1.0, lambda=10 
## + Fold5: alpha=0.0, lambda=10 
## - Fold5: alpha=0.0, lambda=10 
## + Fold5: alpha=0.2, lambda=10 
## - Fold5: alpha=0.2, lambda=10 
## + Fold5: alpha=0.4, lambda=10 
## - Fold5: alpha=0.4, lambda=10 
## + Fold5: alpha=0.6, lambda=10 
## - Fold5: alpha=0.6, lambda=10 
## + Fold5: alpha=0.8, lambda=10 
## - Fold5: alpha=0.8, lambda=10 
## + Fold5: alpha=1.0, lambda=10 
## - Fold5: alpha=1.0, lambda=10

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.

## Aggregating results
## Selecting tuning parameters
## Fitting alpha = 0.2, lambda = 0.1 on full training set

elnet

## glmnet 
## 
## 3794 samples
## 9000 predictors
## 
## No pre-processing
## Resampling: Bootstrapped (5 reps) 
## Summary of sample sizes: 759, 759, 758, 759, 759 
## Resampling results across tuning parameters:
## 
##   alpha  lambda  RMSE       Rsquared   MAE      
##   0.0    1e-03   0.3817416  0.7719830  0.3058639
##   0.0    1e-02   0.3817416  0.7719830  0.3058639
##   0.0    1e-01   0.3817416  0.7719830  0.3058639
##   0.0    1e+00   0.3817416  0.7719830  0.3058639
##   0.0    1e+01   0.3929088  0.7695436  0.3182781
##   0.2    1e-03   0.4144098  0.7191260  0.3157629
##   0.2    1e-02   0.4144098  0.7191260  0.3157629
##   0.2    1e-01   0.3714836  0.7719132  0.2890308
##   0.2    1e+00   0.5455076  0.6625650  0.4716594
##   0.2    1e+01   0.7749574        NaN  0.6959006
##   0.4    1e-03   0.4193101  0.7136845  0.3187730
##   0.4    1e-02   0.4193101  0.7136845  0.3187730
##   0.4    1e-01   0.3767381  0.7718409  0.2988981
##   0.4    1e+00   0.6927642  0.5635876  0.6180105
##   0.4    1e+01   0.7749574        NaN  0.6959006
##   0.6    1e-03   0.4214262  0.7113153  0.3200033
##   0.6    1e-02   0.4184079  0.7147409  0.3178406
##   0.6    1e-01   0.3971945  0.7528825  0.3187922
##   0.6    1e+00   0.7749574        NaN  0.6959006
##   0.6    1e+01   0.7749574        NaN  0.6959006
##   0.8    1e-03   0.4233367  0.7091882  0.3208346
##   0.8    1e-02   0.4118057  0.7222663  0.3123569
##   0.8    1e-01   0.4172564  0.7338352  0.3384660
##   0.8    1e+00   0.7749574        NaN  0.6959006
##   0.8    1e+01   0.7749574        NaN  0.6959006
##   1.0    1e-03   0.4290790  0.7028063  0.3221763
##   1.0    1e-02   0.4085539  0.7261220  0.3079923
##   1.0    1e-01   0.4349320  0.7181826  0.3563351
##   1.0    1e+00   0.7749574        NaN  0.6959006
##   1.0    1e+01   0.7749574        NaN  0.6959006
## 
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0.2 and lambda = 0.1.

Zatem na podstawie naszej 1 próby cross-validacyjnej ( w praktyce warto zrobić więcej, ale ze względu na czytelność 4b) ograniczmy się do 1 ) powinniśmy wybrać α = 0.2 i λ = 0.1.

Rzućmy teraz okiem na błędy na poszczególnych zbiorach walidacji krzyżowej.

elnet$resample

##        RMSE  Rsquared       MAE Resample
## 1 0.3682216 0.7766761 0.2898645    Fold1
## 2 0.3917260 0.7451115 0.2947760    Fold2
## 3 0.3630105 0.7817005 0.2826452    Fold4
## 4 0.3642935 0.7837571 0.2877927    Fold3
## 5 0.3701668 0.7723209 0.2900754    Fold5

Stąd błąd walidacyjny modelu to

mean(elnet$resample[,'RMSE'])

## [1] 0.3714836

A błąd treningowy to

y_train.pred <- predict(elnet, X_train_scaled) # dokonujemy predykcji na zbiorze treningowym
train.error <- sqrt(mean(( unlist(as.vector(y_train)) - y_train.pred)^2))

train.error

## [1] 0.3275472

To wygląda sensownie. Zazwyczaj błąd treningowy jest mniejszy niż błąd zbioru walidacyjnego.

Zadanie 3

grid <- expand.grid(
  mtry = c(200,250,300,350,400),
  ntree = c(2, 3, 4),
  nodesize= c(5, 10, 15, 20))

Narazie nie będziemy wgłębiać się w preprocessing, ze względu na to, że random forest jest metodą “off-the-shelf”, czyli tak naprawdę sam w sobie jest gotowy w użyciu.

set.seed(1)

X_train.ext <- X_train
X_train.ext$y <- y_train

# Bierzemy te same indeksy co przy elasticnet

mse.list <- list() # tu będziemy wrzucać wyniki dla zbiory walidacyjnego dla danego zestawu hiperparametrów


for (i in 1:nrow(grid)){
  
  mse.oneforest <- list() # lista mse dla zbiorów z cv dla bieżącego drzewa
  
  for (j in 1:5){
    

    
    my.randomForest <- randomForest(x=X_train.ext[unlist(flds[j]),c(1:9000)], y= X_train.ext[unlist(flds[j]),c(9001)],
                                    mtry=grid$mtry[i], ntree=grid$ntree[i], nodesize=grid$nodesize[i])
  mse.oneforest <- append(mse.oneforest, tail(my.randomForest$mse, n=1)) # https://stat.ethz.ch/pipermail/r-help/2004-April/049943.html ,                                                          ostatni zwracany mse to mse całego lasu
  }
  mse.list <- append(mse.list, mean(unlist((mse.oneforest)))) # dodajemy średnią z listy dla bieżącego drzewa

}

mse.min.ind <- which.min(mse.list)
grid[mse.min.ind,]

##    mtry ntree nodesize
## 60  400     4       20

Ponieważ w pętli nie przechwytujemy najlepszego modelu to wytrenujmy random forest ponownie dla optymalnych hiperparametrów ( opisywany przykład jest na tyle niewielki, że możemy sobie pozwolić na takie podejście).

final.forest.rmse.cv <- list()

for (j in 1:5){
    
    final.randomForest <- randomForest(x=X_train.ext[unlist(flds[j]),c(1:9000)], y= X_train.ext[unlist(flds[j]),c(9001)],
                                    mtry=grid$mtry[mse.min.ind], ntree=grid$ntree[mse.min.ind], nodesize=grid$nodesize[mse.min.ind])
  final.forest.rmse.cv <- append(final.forest.rmse.cv, sqrt(tail(final.randomForest$mse, n=1))) # https://stat.ethz.ch/pipermail/r-help/2004-April/049943.html ,                                                          ostatni zwracany mse to mse całego lasu, sqrt bo chcemy rmse
  
}

final.forest.rmse.cv

## [[1]]
## [1] 0.5362993
## 
## [[2]]
## [1] 0.4871534
## 
## [[3]]
## [1] 0.4570505
## 
## [[4]]
## [1] 0.5024537
## 
## [[5]]
## [1] 0.5018022

unlist(final.forest.rmse.cv)

## [1] 0.5362993 0.4871534 0.4570505 0.5024537 0.5018022

Tworzymy podsumowanie tabelaryczne

ElasticNet <- unlist(elnet$resample['RMSE'])
   
RandomForest <- unlist(final.forest.rmse.cv)
    
summary.df<- data.frame(ElasticNet, RandomForest)
rownames(summary.df) <- c("Fold1", "Fold2", "Fold3", "Fold4", "Fold5")
summary.df

##       ElasticNet RandomForest
## Fold1  0.3682216    0.5362993
## Fold2  0.3917260    0.4871534
## Fold3  0.3630105    0.4570505
## Fold4  0.3642935    0.5024537
## Fold5  0.3701668    0.5018022

Widzimy, że ElasticNet dla naszych foldów radzi sobie lepiej.

Możemy jeszcze sprawdzić mse na zbiorze treningowym. Dla elasticnet błąd treningowy przechowaliśmy w zmiennej train.error. Teraz policzymy dla RandomForest

elnet.train.error <- train.error

final.randomForest <- randomForest(x=X_train.ext[,c(1:9000)], y= X_train.ext[,c(9001)],
                                    mtry=grid$mtry[mse.min.ind], ntree=grid$ntree[mse.min.ind], nodesize=grid$nodesize[mse.min.ind])

random_forest.train.error <- sqrt(tail(final.randomForest$mse, n=1))

elnet.train.error

## [1] 0.3275472

random_forest.train.error

## [1] 0.4433093

Ponownie ElasticNet radzi sobie na zbiorze lepiej. Zauważmy też, że w poprzednich rozważaniach wybraliśmy model o największych hiperparametrach spośród siatki. Niestety nie mamy dostępu do zbioru testowego. Ale na podstawie tego jak modela radzą sobie na zbiorze walidacyjnym, stwierdzamy, że ElasticNet jest lepszym wyborem, bo błąd jest mniejszy.

Porównajmy jeszcze nasze modele z modelem referencyjnym

reference.rmse.cv <- list()

for (j in 1:5){
    
  avg.y_train <- mean(X_train.ext[unlist(flds[j]),c(9001)])
    
    reference.rmse.cv <- append(reference.rmse.cv, sqrt(mean((X_train.ext[unlist(flds[j]),c(9001)] - avg.y_train)^2)))
  
}

Reference <- unlist(reference.rmse.cv)

summary.df$ReferenceModel <- Reference

summary.df

##       ElasticNet RandomForest ReferenceModel
## Fold1  0.3682216    0.5362993      0.7783484
## Fold2  0.3917260    0.4871534      0.7762131
## Fold3  0.3630105    0.4570505      0.7669997
## Fold4  0.3642935    0.5024537      0.7743854
## Fold5  0.3701668    0.5018022      0.7782463

Jak widzimy model referencyjny wypada słabo i tego mogliśmy się spodziewać. Ale warto było kontrolnie sprawdzić jego wyniki. Gdyby jego wyniki były lepsze od któregoś z modeli, to prawie na pewno popełnilibyśmy po drodze jakiś błąd.

Zadanie 4

W tym zadaniu podejścia mogą być różne. W szczególności warto zastanowić się nad tym jakie modele znamy i czy możemy je wykorzystać. W poprzednich podpunktach rzuciliśmy okiem na elasticnet i randomforest. Patrząc na wyniki możemy stwierdzić, że randomforest nie poradził sobie najlepiej, jak na złożoność jego modelu. Wyświetlmy informacje o wytrenowanym modelu:

print(final.randomForest)

## 
## Call:
##  randomForest(x = X_train.ext[, c(1:9000)], y = X_train.ext[,      c(9001)], ntree = grid$ntree[mse.min.ind], mtry = grid$mtry[mse.min.ind],      nodesize = grid$nodesize[mse.min.ind]) 
##                Type of random forest: regression
##                      Number of trees: 4
## No. of variables tried at each split: 400
## 
##           Mean of squared residuals: 0.1965231
##                     % Var explained: 67.27

Naszą uwagę powinien przykuć dość niski procent wyjaśnionych danych. Wiemy, że czas naszej pracy jest skończony, dlatego musimy obrać taktykę. Zdecyduje się na trenowanie randomforest na lepszym gridzie, bo spodziewam się po tym modelu, że może poradzić sobie znacznie lepiej. Jeśli otrzymam zadawalające wyniki to zakończę rozważania, jeśli nie, to zmienię strategię.

Powtarzamy nasze rozważania dla random forest, ale dla siatki o większych parametrach.

grid <- expand.grid(
  mtry = c(1000,5000),
  ntree = c(500,750))

set.seed(1)

mse.list <- list() # tu będziemy wrzucać wyniki dla zbiory walidacyjnego dla danego zestawu hiperparametrów


for (i in 1:nrow(grid)){
  
  mse.oneforest <- list() # lista mse dla zbiorów z cv dla bieżącego drzewa
  
  for (j in 1:5){
    

    
    my.randomForest <- randomForest(x=X_train.ext[unlist(flds[j]),c(1:9000)], y= X_train.ext[unlist(flds[j]),c(9001)],
                                    mtry=grid$mtry[i], ntree=grid$ntree[i], nodesize=grid$nodesize[i])
  mse.oneforest <- append(mse.oneforest, tail(my.randomForest$mse, n=1)) # https://stat.ethz.ch/pipermail/r-help/2004-April/049943.html ,                                                          ostatni zwracany mse to mse całego lasu
  }
  mse.list <- append(mse.list, mean(unlist((mse.oneforest)))) # dodajemy średnią z listy dla bieżącego drzewa

}

mse.min.ind <- which.min(mse.list)
grid[mse.min.ind,]

##   mtry ntree
## 2 5000   500

Zatem wytrenujmy model dla uzyskanych hiperparametrów, na całych danych.

( Kod gdyby google colab nie chcial współpracować )

# set.seed(1)
# 
# final.randomForest <- randomForest(x=X_train, y=y_train, mtry=5000, ntree=500)
# final.forest.rmse.cv <- sqrt(tail(final.randomForest$mse, n=1)) # https://stat.ethz.ch/pipermail/r-help/2004-April/049943.html , ostatni zwracany mse to mse całego lasu, sqrt bo chcemy rmse
# 
# print(final.forest.rmse.cv)
# print(final.randomForest)

Z powodu słabych podzespołów mojego komputera, skorzystałem z google colab:

https://colab.research.google.com/drive/1RART58R4S0Nn-FYAMZDNJXmK7TnJXrEy?usp=sharing

Model trenował się +-4h, uzyskane wyniki są zadawalające ( 83% wyjaśnionych danych i całkiem dobry bląd dla danych testowych). Potencjalnie moglibyśmy się jeszcz pokusić o spróbowanie grida o większych wartościach mtry.

( Kod do predykcji gdyby google colab nie chcial współpracować )

# final.predictions <- predict(final.randomForest, newdata = X_test)
# Id <- c(0:669)
# final.df <- data.frame(Id,final.predictions)
# 
# colnames(final.df) <- c('Id','Expected')
#                   
# write.csv(final.df, file="predictions.csv", row.names = FALSE)

Bibliografia: - Wykłady i laboratoria przedmiotu Statystyczna Analiza Danych wydziału MIMUW: https://usosweb.mimuw.edu.pl/kontroler.php?_action=katalog2/przedmioty/pokazPrzedmiot&kod=1000-714SAD

• https://stats.stackexchange.com/questions/132652/how-to-determine-which-distribution-fits-my-data-best
• https://stackoverflow.com/questions/44507568/error-when-fitting-a-beta-distribution-the-function-mle-failed-to-estimate-the
• https://stats.stackexchange.com/questions/126539/testing-whether-data-follows-t-distribution/126552#126552
• https://www.rdocumentation.org/packages/Sim.DiffProc/versions/2.5/topics/test_ks_dbeta
• https://stackoverflow.com/questions/2661402/given-a-set-of-random-numbers-drawn-from-a-continuous-univariate-distribution-f
• https://towardsdatascience.com/selecting-the-best-predictors-for-linear-regression-in-r-f385bf3d93e9
• https://gdudek.el.pcz.pl/files/SUS/SUS_wyklad6.pdf
• https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html
• https://machinelearningmastery.com/elastic-net-regression-in-python/
• https://www.datacamp.com/tutorial/tutorial-ridge-lasso-elastic-net
• https://towardsdatascience.com/optimization-of-elastic-net-regularization-for-predicting-soil-clay-content-7df09f6a68f8
• https://corporatefinanceinstitute.com/resources/knowledge/other/elastic-net/
• https://daviddalpiaz.github.io/r4sl/elastic-net.html
• https://datascience.stackexchange.com/questions/28158/how-to-calculate-the-fold-number-k-fold-in-cross-validation
• https://www.baeldung.com/cs/train-test-datasets-ratio
• https://www.analyticsvidhya.com/blog/2022/02/k-fold-cross-validation-technique-and-its-essentials/
• https://scikit-learn.org/stable/modules/cross_validation.html
• https://www.projectpro.io/recipes/create-and-optimize-baseline-elasticnet-regression-model
• https://www.youtube.com/watch?v=1dKRdX9bfIo&ab_channel=StatQuestwithJoshStarmer
• https://stats.stackexchange.com/questions/187335/validation-error-less-than-training-error
• https://rpubs.com/phamdinhkhanh/389752
• https://stats.stackexchange.com/questions/172842/best-practices-with-data-wrangling-before-running-random-forest-predictions